# This file estimates longterm effects of twinning for mothers who are 20 years older 

# install.packages("dplyr")
# install.packages("date")
# install.packages("AER")

library(dplyr)
library(date)
library(AER)

# Load data
load("mothers_old_20.rdata")

# Run OLS models for any children 

ols_dat1 <- 
  mothers %>% 
  mutate(children = ! is.na(age_oldest)) %>%
  filter(! is.na(stemt)) %>%
  mutate(vote = stemt == "Stemte")

ols_dat1 <-
  within(ols_dat1, {
    age        <- ntile(date, 25)
  })

ols_mod1 <- lm(vote ~ children + as.factor(age), data = ols_dat1)

# Filter of if age of oldest is unkwown (select parents )
mothers <-
  mothers %>%
  filter(!is.na(age_oldest)) 

# split data file to each election 
# 2009 election
# Subset and create vote variable
# filter of parents with children born after the election

mothers09 <- 
  mothers %>%
  filter(! is.na(stemte_2009) & first09 > 0) %>%
  filter(age_oldest <= (14565 - 20*365 - 5)) %>%
  mutate(vote = stemte_2009 == "Stemte",
         no_children = no_children_09)

#create age dummies for mothers 
#dichotomize twin first
mothers09 <-
  within(mothers09, {
    age        <- ntile(date, 25)
    twin_first <- twin_first > 0
  })

# create dummies for age of children grouped by age of mother dummies
mothers09 <-
  mothers09 %>%
  group_by(age) %>%
  mutate(child_age = ntile(age_oldest, 4))

## Repeat for 13
mothers13 <- 
  mothers %>%
  filter(! is.na(stemt) & first13 > 0) %>%
  filter(age_oldest <= 16028 - 20*365 - 5) %>%
  mutate(vote = stemt == "Stemte",
         no_children = no_children_13)

mothers13 <-
  within(mothers13, {
    age        <- ntile(date, 25)
    twin_first <- twin_first > 0
  })

mothers13 <-
  mothers13 %>%
  group_by(age) %>%
  mutate(child_age = ntile(age_oldest, 4))

## Repeat for 14
mothers14 <- 
  mothers %>%
  filter(! is.na(ep_stemt) & first13 > 0) %>%
  filter(age_oldest <= 16215 - 20*365 - 5) %>%
  mutate(vote = ep_stemt,
         no_children = no_children_14)

mothers14 <-
  within(mothers14, {
    age        <- ntile(date, 25)
    twin_first <- twin_first > 0
  })

mothers14 <-
  mothers14 %>%
  group_by(age) %>%
  mutate(child_age = ntile(age_oldest, 4))

# Run models 
# Simple models. Biased but used for turnout of "control" group
mod09_r <- lm(vote ~ twin_first, data = mothers09)
mod13_r <- lm(vote ~ twin_first, data = mothers13)
mod14_r <- lm(vote ~ twin_first, data = mothers14)

# Models with fixed effect for mothers and children's age
mod09 <- lm(vote ~ twin_first + as.factor(age) * as.factor(child_age), data = mothers09)
mod13 <- lm(vote ~ twin_first + as.factor(age) * as.factor(child_age), data = mothers13)
mod14 <- lm(vote ~ twin_first + as.factor(age) * as.factor(child_age), data = mothers14)


## Instrument number of children with first child being twin
mod09_iv  <- ivreg(vote ~ no_children + as.factor(age) * as.factor(child_age) | 
                     twin_first + as.factor(age) * as.factor(child_age), data = mothers09)
mod13_iv  <- ivreg(vote ~ no_children + as.factor(age) * as.factor(child_age) | 
                     twin_first + as.factor(age) * as.factor(child_age), data = mothers13)
mod14_iv  <- ivreg(vote ~ no_children + as.factor(age) * as.factor(child_age) | 
                     twin_first + as.factor(age) * as.factor(child_age), data = mothers14)

## Run first-stage models 
mod09_fs <- lm(no_children ~ twin_first + as.factor(age) * as.factor(child_age), data = mothers09)
mod13_fs <- lm(no_children ~ twin_first + as.factor(age) * as.factor(child_age), data = mothers13)
mod14_fs <- lm(no_children ~ twin_first + as.factor(age) * as.factor(child_age), data = mothers14)

## Run simple OLS for relationship between no of children and turnout controlling for age fixed effects 

mothers13 <-
  within(mothers13, {
    no_children_ind <- no_children
    no_children_ind[no_children_ind > 5] <- 5
  })

ols_mod2 <- 
  lm(vote ~ as.factor(no_children_ind) + as.factor(age) * as.factor(child_age), data = mothers13)

# Compile results 

table_ols_13 <- 
  rbind(summary(ols_mod1)$coefficients[1:2,],
        summary(ols_mod2)$coefficients[2:5,])

table_mom_09 <-
  rbind(summary(mod09_r)$coefficients,
        summary(mod09)$coefficients[2,],
        summary(mod09_iv)$coefficients[2,])

table_mom_13 <-
  rbind(summary(mod13_r)$coefficients,
        summary(mod13)$coefficients[2,],
        summary(mod13_iv)$coefficients[2,])

table_mom_14 <-
  rbind(summary(mod14_r)$coefficients,
        summary(mod14)$coefficients[2,],
        summary(mod14_iv)$coefficients[2,])

ns <- rbind(
  with(mothers09, table(twin_first)),
  with(mothers13, table(twin_first)),
  with(mothers14, table(twin_first)),
  with(ols_dat1, table(children))
)

tab_first_stage <- 
  rbind(summary(mod09_fs)$coefficients[2,],
        summary(mod13_fs)$coefficients[2,],
        summary(mod14_fs)$coefficients[2,])

outtab_mom <- 
  list(table_ols_13,
       table_mom_09,
       table_mom_13,
       table_mom_14,
       ns,
       with(mothers13, table(no_children_ind)),
       tab_first_stage)


save(outtab_mom, file = "mother_twin_old_20_res.rdata")

